使用Python绘制原子轨道与杂化轨道(的等概率密度面), 效果如下:
此处演示使用的版本为Python 3.9.7, 建议版本号>= 3.7.
其他还会使用的包及版本号如下: (版本号有一点差别也不要紧, 通常不会有大问题)
numpy 1.21.4scipy 1.10.1scikit-image 0.21.0matplotlib 3.5.0其中$\hbar$为约化普朗克常数, $\mu={m_\mathrm{p} m_\mathrm{e} \over m_\mathrm{p} + m_\mathrm{e}}$为原子核与电子的约化质量, $\nabla^2$为拉普拉斯算符(梯度的散度), $e$为电子电量, $\varepsilon_0$为真空介电常数, $E$为原子轨道能量.
为了让$\Psi$有物理意义(单值、有限), 求解过程还会自然地引入三个量子数: 主量子数$n$, 角量子数$\ell$和磁量子数$m$. 它们需要满足的关系有:
$n \in \mathbb{N}^{+},\, \ell \in \mathbb{N},\, m \in \mathbb{Z},\, n > \ell,\,\ell \geq |m|$.
不同的$(n, \ell, m)$组合产生不同的原子轨道, 如$(n, \ell, m) = (1,0,0)$表示$\text{1s}$轨道, $(n, \ell, m) = (3, 2, 1)$表示$\text{3d}_{1}$轨道.
其中$\rho = {2 r \over n a_{0}^{*}}$, 是一个简化书写形式的无量纲中间量. $a_{0}^{*}$是约化玻尔半径, 大约是$52.9 \,\mathrm{pm}$. $L_{n-\ell-1}^{2\ell + 1} \left( \cdot \right)$是连带拉盖尔多项式. $Y_{\ell}^{m}(\theta, \varphi)$是球谐函数.
为了方便画图, 这里采用类似原子单位的处理方式, 令$a_{0}^{*} = 1$. 此时$r$代表着真实的径向距离是$a_{0}^{*} (= 52.9 \,\mathrm{pm})$的多少倍.
注: 此处使用的球坐标系的形式——径向距离$r$, 极角$\theta$和方位角$\varphi$.
需要提前安装numpy, scipy, scikit-image, matplotlib.
导入的主要内容有:
marching cubes算法import numpy as np
# 特殊函数:阶乘(factorial)、连带拉盖尔多项式(genlaguerre)、球谐函数(sph_harm)
from scipy.special import factorial, genlaguerre, sph_harm
# 等值面搜索:marching_cubes
from skimage.measure import marching_cubes
# 绘图组件
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
# Jupyter Notebook内嵌绘图用,本地运行则不写这一行
%matplotlib notebook
把定义抄过来:
\begin{equation} \Psi_{n \ell m} (r, \theta, \varphi) = \sqrt{\left( {2 \over n a_{0}^{*}} \right)^3 {(n - \ell - 1)! \over 2 n (n+\ell)!}} {\rho}^{\ell} \exp\left( -{\rho \over 2} \right) L_{n-\ell-1}^{2\ell + 1} \left( \rho \right) \cdot Y_{\ell}^{m}(\theta, \varphi) \end{equation}其中$\rho = {2 r \over n a_{0}^{*}}$, 是一个简化书写形式的无量纲中间量. $a_{0}^{*}$是约化玻尔半径, 大约是$52.9 \,\mathrm{pm}$. $L_{n-\ell-1}^{2\ell + 1} \left( \cdot \right)$是连带拉盖尔多项式. $Y_{\ell}^{m}(\theta, \varphi)$是球谐函数.
为了方便画图, 这里采用类似原子单位的处理方式, 令$a_{0}^{*} = 1$. 此时$r$代表着真实的径向距离是$a_{0}^{*} (= 52.9 \,\mathrm{pm})$的多少倍.
注: 此处使用的球坐标系的形式——径向距离$r$, 极角$\theta$和方位角$\varphi$.
直接翻译成代码就行.
# 定义约化玻尔半径a0, 原子单位下直接设置为1
a0 = 1
def hydrogen_wave_function(n, l, m):
# 径向部分R(r), 注意genlaguerre的参数顺序
def R(r):
factor = np.sqrt((2. / (n * a0))**3 * \
factorial(n - l - 1) / (2 * n * factorial(n + l))) # 系数
rho = 2 * r / (n * a0) # 中间量rho
return factor * (rho**l) * np.exp(-rho / 2) * \
genlaguerre(n - l - 1, 2 * l + 1)(rho)
# 角向部分Y(theta, phi)就是球谐函数sph_harm
# 注意量子数l与m 和 自变量theta与phi的顺序都是颠倒的!
def Y(theta, phi):
return sph_harm(m, l, phi, theta)
# 径向部分R(r)和角向部分Y(theta, phi)相乘,得到一个关于(r, theta, phi)的函数,即原子轨道
return lambda r, theta, phi: R(r) * Y(theta, phi)
再代入$n$, $\ell$和$m$三个量子数就得到了对应的轨道波函数hydrogen_wave_function(n, l, m).
以$\text{2p}_0$轨道(又名$\text{2p}_z$轨道)为例, 对应的量子数为$(n, \ell, m) = (2, 1, 0)$:
# 原子轨道参数,当前为2p0轨道
n, l, m = 2, 1, 0
# 2p0轨道波函数
psi = hydrogen_wave_function(n, l, m)
# `psi`本身是一个函数,后面还可以继续带参数(r, theta, phi)来计算实际的波函数值
print(f"{psi = }")
print(f"{psi(1, np.pi/4, np.pi/3) = }") # 注意结果是复数
psi = <function hydrogen_wave_function.<locals>.<lambda> at 0x0000020F8F9500D0> psi(1, np.pi/4, np.pi/3) = (0.04277478503902706+0j)
把空间网格化,确定边界limit和三个正方向上取的点数n_points.
网格化的区域是一个边长为2 * limit且中心在原点的正方体. 每个方向上均匀取n_points个点, 网格点总数为(n_points)^3.
# limit: 确定网格点的区间范围
limit = 10
# n_points: 每个方向上均匀取点的数目
n_points = 50
# vec: 临时变量, 从 -limit 到 limit 均匀取 n_points 个点得到的向量
vec = np.linspace(-limit, limit, n_points)
print(f"{vec = }")
print(f"{len(vec) = }")
vec = array([-10. , -9.59183673, -9.18367347, -8.7755102 ,
-8.36734694, -7.95918367, -7.55102041, -7.14285714,
-6.73469388, -6.32653061, -5.91836735, -5.51020408,
-5.10204082, -4.69387755, -4.28571429, -3.87755102,
-3.46938776, -3.06122449, -2.65306122, -2.24489796,
-1.83673469, -1.42857143, -1.02040816, -0.6122449 ,
-0.20408163, 0.20408163, 0.6122449 , 1.02040816,
1.42857143, 1.83673469, 2.24489796, 2.65306122,
3.06122449, 3.46938776, 3.87755102, 4.28571429,
4.69387755, 5.10204082, 5.51020408, 5.91836735,
6.32653061, 6.73469388, 7.14285714, 7.55102041,
7.95918367, 8.36734694, 8.7755102 , 9.18367347,
9.59183673, 10. ])
len(vec) = 50
得到单个方向上的均匀取点后, 用np.meshgrid生成三维网格的格点坐标.
# 生成三维网格
# X记录每个点的x坐标,Y记录y坐标, Z记录z坐标
X, Y, Z = np.meshgrid(vec, vec, vec)
print(f"{X.shape = }", f"{Y.shape = }", f"{Z.shape = }")
X.shape = (50, 50, 50) Y.shape = (50, 50, 50) Z.shape = (50, 50, 50)
# 查看每个点的直角坐标(x, y, z)
# 后续并不使用`coord_xyz`,只用来确认数据
coords_xyz = np.vstack(list(map(np.ravel, (X, Y, Z)))).T
print(f"{len(coords_xyz) = }")
print(coords_xyz)
len(coords_xyz) = 125000 [[-10. -10. -10. ] [-10. -10. -9.59183673] [-10. -10. -9.18367347] ... [ 10. 10. 9.18367347] [ 10. 10. 9.59183673] [ 10. 10. 10. ]]
或者直接把它们标出来:
格点是直角坐标$(x,\, y,\, z)$, 然而$\Psi$是关于$(r,\, \theta,\, \varphi)$的函数, 需要把格点按下式变换为球坐标.
\begin{equation} \left\{ \begin{aligned} x & = r \sin\theta \cos\varphi \\\\ y & = r \sin\theta \sin\varphi \\\\ z & = r \cos\theta \end{aligned} \right. \quad \Leftrightarrow \quad \left\{ \begin{aligned} r & = \sqrt{x^2 + y^2 + z^2} \\ \theta & = \arccos {z \over r} = \arccos {z \over \sqrt{x^2 + y^2 + z^2}} \\ \varphi & = \arctan {y \over x}\\ \end{aligned} \right. \end{equation}# 由坐标变换,每个格点从直角坐标系下的(X, Y, Z)变成球坐标系下的(R, THETA, PHI)
R = np.sqrt(X**2 + Y**2 + Z**2)
THETA = np.arccos(Z / R)
PHI = np.arctan2(Y, X)
print(f"{R.shape = }", f"{THETA.shape = }", f"{PHI.shape = }")
R.shape = (50, 50, 50) THETA.shape = (50, 50, 50) PHI.shape = (50, 50, 50)
# 查看每个点的球坐标(r, theta, phi)
# 后续并不使用`coord_rtf`,只用来确认数据
coords_rtf = np.vstack(list(map(np.ravel, (R, THETA, PHI)))).T
print(f"{len(coords_rtf) = }")
print(coords_rtf)
len(coords_rtf) = 125000 [[17.32050808 2.18627604 -2.35619449] [17.08810498 2.16677211 -2.35619449] [16.86237997 2.14673823 -2.35619449] ... [16.86237997 0.99485442 0.78539816] [17.08810498 0.97482054 0.78539816] [17.32050808 0.95531662 0.78539816]]
# 代入波函数`psi`中, 得到每个格点的波函数值
psi_values = psi(R, THETA, PHI)
print(f"{psi_values.shape = }")
print(f"{psi_values[0][0][0] = }") # 直角坐标系下(-10, -10, -10)处的波函数值
psi_values.shape = (50, 50, 50) psi_values[0][0][0] = (-0.0001728819028298787+0j)
使用marching_cubes寻找概率密度$|\Psi|^2$的等值面(i.e. 所有$|\Psi|^2=\text{常数}C$的点的集合)
# prob_dens: 概率密度|Psi|^2
prob_dens = np.abs(psi_values)**2
# iso_value: 目标值,也就是上面的常数C
iso_value = 4e-4
# 初次尝试
verts, faces, _, _ = marching_cubes(
prob_dens,
level=iso_value,
)
print(f"{verts = }")
print(f"{faces = }")
verts = array([[13.979592, 21. , 16. ],
[14. , 21. , 15.840323],
[14. , 20.949831, 16. ],
...,
[35.204258, 27. , 34. ],
[35.02041 , 28. , 16. ],
[35.02041 , 28. , 33. ]], dtype=float32)
faces = array([[ 2, 1, 0],
[ 0, 3, 2],
[ 6, 5, 4],
...,
[3830, 3902, 3813],
[3903, 3836, 3822],
[3836, 3903, 3825]])
结果为一个三角化的等值面. verts本应该代表顶点坐标, 但是数值不正确(默认从0开始计数, 且步长默认为1). faces代表每个三角面的顶点编号.
为了获得正确的顶点坐标, 需要再做三件事:
marching_cubes加上spacing参数表示每个方向的步长. 步长step可以根据limit和n_points参数简单计算得到.[-limit, limit]是对称的, 所以直接给每个坐标减掉向量(limit, limit, limit)即可.marching_cubes坐标输出顺序的影响, 所有顶点的$x$坐标和$y$坐标需要互换(相关说明见附录A, 这里只提结论, 不再展开)# 计算步长step: [-limit, limit]区间被分成(n_points-1)份,每段的长度就是step
step = 2 * limit / (n_points - 1)
# 把step传入marching_cubes的spacing参数中
verts, faces, _, _ = marching_cubes(
prob_dens,
level=iso_value,
spacing=(step, step, step),
)
# 将结果平移至原点
verts -= limit
# 最后互换x坐标和y坐标
verts[:, [0, 1]] = verts[:, [1, 0]]
# 结果确认
print(f"{verts = }")
verts = array([[-1.42857143, -4.29404395, -3.46938776],
[-1.42857143, -4.28571429, -3.53456186],
[-1.44904857, -4.28571429, -3.46938776],
...,
[ 1.02040816, 4.36908488, 3.87755102],
[ 1.42857143, 4.29404434, -3.46938776],
[ 1.42857143, 4.29404434, 3.46938776]])
此时的verts坐标全部正常.
获得正确的等值面后, 接下来就是把它画出来. 这里介绍3种可视化方法:
ax.plot_trisurfPoly3DCollectionobj文件中, 借助其他软件(如Blender, 3D Viewer等)展示结果为了前两种方法的方便, 先写个辅助函数, 调用时返回一个调整好的画布和内框.
# 辅助函数,每次返回一个调整好的画布`fig`和内框`ax`
def new_fig_and_ax(plot_range=10):
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection="3d")
# 设置坐标比例关系,让图像变得方正
ax.set_box_aspect([1, 1, 1])
# 设置合适的视图俯仰角和方位角,让坐标轴方向符合直观
# x轴对着屏幕外方向,y轴朝右,z轴朝上
ax.view_init(elev=30, azim=20)
# 设置坐标轴标签
ax.set_xlabel("$x$")
ax.set_ylabel("$y$")
ax.set_zlabel("$z$")
# 设置坐标轴范围为[-plot_range, plot_range]
ax.set_xlim(-plot_range, plot_range)
ax.set_ylim(-plot_range, plot_range)
ax.set_zlim(-plot_range, plot_range)
# 设置坐标轴刻度
ticks = np.linspace(-plot_range, plot_range, 5)
ax.set_xticks(ticks)
ax.set_yticks(ticks)
ax.set_zticks(ticks)
return fig, ax
ax.plot_trisurf¶此处plot_trisurf的传参顺序: 所有点的X坐标, 所有点的Y坐标, 所有三角面的顶点编号, 所有点的Z坐标.
上述参数3来自faces; 其他参数来自verts, 但需要把同类坐标整合在一起.
# 方法一:ax.plot_trisurf
fig, ax = new_fig_and_ax()
iso_surface = ax.plot_trisurf(
verts[:, 0], verts[:, 1], faces, verts[:, 2],
lw=0, # 线宽设置为0
cmap="coolwarm_r", # 可以带上颜色让图像更好看,但是注意这个颜色与波函数的物理意义无关
)
# 展示结果
plt.show()
# 或者使用plt.savefig(...)把图像存储到文件中
# plt.savefig(
# "orbital_plot.png",
# transparent=True,
# bbox_inches='tight',
# dpi=600,
# )
Poly3DCollection¶结果和plot_trisurf类似, 但是Poly3DCollection的参数更简单: 每个三角面三个顶点坐标的集合, 只需要用verts[faces]即可.
# 方法二:Poly3DCollection
fig, ax = new_fig_and_ax()
mesh = Poly3DCollection(verts[faces], lw=0.1) # 线宽0.1
mesh.set_facecolor("yellow") # 黄色的表面
mesh.set_edgecolor("grey") # 灰色的边线
ax.add_collection3d(mesh) # 直接加入ax中
plt.show()
obj文件, 借助其他软件展示结果¶处理起来最自由的一种方式, 而且obj文件的结构非常简单, 这里只需要知道:
o开头的行 —— 对象信息: o 【对象名称】v开头的行 —— 顶点信息: v 【顶点x坐标】 【顶点y坐标】 【顶点z坐标】f开头的行 —— 面信息: f 【三角面顶点1编号】【三角面顶点2编号】【三角面顶点3编号】# 方法三: 写入obj文件
filename = "isosurface.obj"
with open(filename, "w") as f:
# (1) o开头的行:写入对象名称,可以自行指定
f.write(f"o Isosurface-2pz {iso_value}\n")
# (2) v开头的行:写入顶点坐标(x, y, z)
for vert_coords in verts:
x, y, z = vert_coords
f.write("v %s %s %s \n" % (x, y, z))
# (3) f开头的行:写入三角面的顶点编号(id_1, id_2, id_3),注意编号从1开始计
for vert_ids in faces:
id_1, id_2, id_3 = vert_ids + 1
f.write("f %s %s %s \n" % (id_1, id_2, id_3))
# obj文件内容
with open(filename, "r") as f:
lines = f.readlines()
print(f"{len(lines) = }\n")
print("".join(lines[:5]), "...")
print("".join(lines[-5:]))
len(lines) = 11705 o Isosurface-2pz 0.0004 v -1.4285714285714288 -4.294043949672154 -3.4693877551020407 v -1.4285714285714288 -4.285714285714286 -3.534561857885244 v -1.4490485677913743 -4.285714285714286 -3.4693877551020407 v -1.4285714285714288 -4.285714285714286 -3.2705976525131533 ... f 3825 3902 3805 f 3903 3831 3811 f 3831 3903 3814 f 3904 3837 3823 f 3837 3904 3826
然后将obj文件导入其他软件中处理, 例如下面是用Blender绘制的结果.
因为波函数的角向部分(也就是球谐函数)含有$\mathrm{e}^{\mathrm{i} m \varphi}$, 其中$\mathrm{i}=\sqrt{-1}$为虚数单位, 所以刚才的波函数大多是复函数. 不过, 我们可以选择$\Psi_{n \ell m}$和$\Psi_{n \ell \bar{m}}$(角标$\bar{m}$代表$-m$)的线性组合来消去复数, 从而得到实函数. 原理就是著名的欧拉公式:
\begin{equation} {1 \over 2} \left( \mathrm{e}^{\mathrm{i} m \varphi} + \mathrm{e}^{-\mathrm{i} m \varphi} \right) = \cos{m \varphi} \\ {1 \over 2\mathrm{i}} \left( \mathrm{e}^{\mathrm{i} m \varphi} - \mathrm{e}^{-\mathrm{i} m \varphi} \right) = \sin{m \varphi} \end{equation}如果$\Psi_{n \ell m}$和$\Psi_{n \ell \bar{m}}$的角向部分恰好互为共轭, 也就是满足$Y_{\ell}^{-m}(\theta, \varphi) = \left[ Y_{\ell}^{m}(\theta, \varphi) \right]^{*}$, 那么就可以直接组合:
\begin{equation} \begin{bmatrix} \Psi_{\text{new}, 1} \\ \Psi_{\text{new}, 2} \end{bmatrix} = {1 \over \sqrt{2}} \begin{bmatrix} 1 & 1 \\ -\mathrm{i} & \mathrm{i} \end{bmatrix} \begin{bmatrix} \Psi_{n \ell m} \\ \Psi_{n \ell \bar{m}} \end{bmatrix} \end{equation}但是scipy使用的球谐函数带有$\text{Condon-Shortley}$相位因子$(-1)^m$, 此时满足的关系则是
$Y_{\ell}^{-m}(\theta, \varphi) = (-1)^m \left[ Y_{\ell}^{m}(\theta, \varphi) \right]^{*}$, 并不能按上式直接加和. 为了消除$\text{CS}$相因子的影响, 当磁量子数$m$为正数时, 我们给它再乘上一个$(-1)^m$. 所以真正的组合其实是这样的:
例如$\text{2p}_{x}$和$\text{2p}_{y}$轨道, 它们是$\Psi_{211}$和$\Psi_{21 \bar{1}}$的线性组合. 将$(n, \ell, m) = (2, 1, 1)$代入上式, 得到:
\begin{equation} \begin{bmatrix} \Psi_{\text{2p}_{x}} \\ \Psi_{\text{2p}_{y}} \end{bmatrix} = {1 \over \sqrt{2}} \begin{bmatrix} -1 & 1 \\ \mathrm{i} & \mathrm{i} \end{bmatrix} \begin{bmatrix} \Psi_{2 1 1} \\ \Psi_{2 1 \bar{1}} \end{bmatrix} \end{equation}# 2px轨道波函数
psi_2px = lambda r, theta, phi: \
(-1) * 1 / np.sqrt(2) * hydrogen_wave_function(2, 1, 1)(r, theta, phi) + \
1 / np.sqrt(2) * hydrogen_wave_function(2, 1, -1)(r, theta, phi)
# 2py轨道波函数
psi_2py = lambda r, theta, phi: \
1j / np.sqrt(2) * hydrogen_wave_function(2, 1, 1)(r, theta, phi) + \
1j / np.sqrt(2) * hydrogen_wave_function(2, 1, -1)(r, theta, phi)
#(思考:这里同时用到了两个波函数,并且计算的结果仍旧是复数,是否有更简洁&更适当的写法?)
#(提示:使用z.real或者z.imag获取复数z的实部或虚部)
按之前的做法画$\text{2p}_{x}$轨道的图像,网格和等值面的常数$C$用之前的值,其他变量名加上2px后缀做区分.
# 找到等值面
iso_value = 4e-4
psi_2px_values = psi_2px(R, THETA, PHI)
prob_dens_2px = np.abs(psi_2px_values)**2 # 2px轨道概率密度
verts_2px, faces_2px, _, _ = marching_cubes(
prob_dens_2px, level=iso_value,
spacing = (step, step, step),
)
verts_2px -= limit
verts_2px[:, [0, 1]] = verts_2px[:, [1, 0]]
# 绘制等值面
fig, ax = new_fig_and_ax()
mesh_2px = Poly3DCollection(verts_2px[faces_2px], lw=0.1) # 线宽0.1
mesh_2px.set_facecolor("yellow") # 黄色的表面
mesh_2px.set_edgecolor("grey") # 灰色的边线
ax.add_collection3d(mesh_2px) # 直接加入ax中
plt.show()
形状与$\text{2p}_{z}$轨道一致, 只是方向沿着$x$轴, 确实是期望的结果.
$\text{2p}_{y}$轨道的画法类似, 不再赘述.
杂化轨道本质上也是原子轨道的线性组合, 所以绘制方法同上面的第二章第5节, 但前提是需要知道组合形式.
以$\text{2s}$, $\text{2p}_{x}$, $\text{2p}_{y}$和$\text{2p}_{z}$轨道构成的$\text{sp}^{3}$杂化轨道为例, 对应的线性组合为:
\begin{equation} \begin{bmatrix} \Psi_{\text{sp}^{3}, a} \\ \Psi_{\text{sp}^{3}, b} \\ \Psi_{\text{sp}^{3}, c} \\ \Psi_{\text{sp}^{3}, d} \end{bmatrix} = {1 \over 2} \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 1 & -1 & -1 \\ 1 & -1 & 1 & -1 \\ 1 & -1 & -1 & 1 \end{bmatrix} \begin{bmatrix} \Psi_{\text{2s}} \\ \Psi_{\text{2p}_{x}} \\ \Psi_{\text{2p}_{y}} \\ \Psi_{\text{2p}_{z}} \end{bmatrix} \end{equation}# 2s轨道
psi_2s = hydrogen_wave_function(2, 0, 0)
# 2px, 2py轨道定义同上,只是这里换一种简洁的写法
psi_2px = lambda r, theta, phi: \
(-1) * np.sqrt(2) * hydrogen_wave_function(2, 1, 1)(r, theta, phi).real
psi_2py = lambda r, theta, phi: \
(-1) * np.sqrt(2) * hydrogen_wave_function(2, 1, 1)(r, theta, phi).imag
# 2pz轨道是最开始介绍的2p0轨道
psi_2pz = hydrogen_wave_function(2, 1, 0)
# 上述杂化轨道sp3,c的定义
psi_sp3_c = lambda r, theta, phi: \
1 / 2 * psi_2s(r, theta, phi) - \
1 / 2 * psi_2px(r, theta, phi) + \
1 / 2 * psi_2py(r, theta, phi) - \
1 / 2 * psi_2pz(r, theta, phi)
# 找到等值面
iso_value = 4e-4
psi_sp3_c_values = psi_sp3_c(R, THETA, PHI)
prob_dens_sp3_c = np.abs(psi_sp3_c_values)**2
verts_sp3_c, faces_sp3_c, _, _ = marching_cubes(
prob_dens_sp3_c, level=iso_value,
spacing = (step, step, step),
)
verts_sp3_c -= limit
verts_sp3_c[:, [0, 1]] = verts_sp3_c[:, [1, 0]]
# 绘制等值面
fig, ax = new_fig_and_ax()
mesh_sp3_c = Poly3DCollection(verts_sp3_c[faces_sp3_c], lw=0.1)
mesh_sp3_c.set_facecolor("cyan")
mesh_sp3_c.set_edgecolor("grey")
ax.add_collection3d(mesh_sp3_c)
ax.view_init(elev=10, azim=50) # 适当调整视角
plt.show()
其余杂化轨道的绘制方法类似, 同样不再赘述.
利用matplotlib(即前两种可视化方法)重新绘制$\text{2p}_{z}$轨道, 但是要使用不同的颜色体现出波函数的正负情况,效果如下:
(提示: 分别获取波函数为正和波函数为负的部分)
$\text{3d}$轨道有两种常见的表现形式, 一种是$(\text{3d}_{-2} ,\, \text{3d}_{-1} ,\, \text{3d}_{0} ,\, \text{3d}_{1} ,\, \text{3d}_{2})$, 另一种是$(\text{3d}_{xy} ,\, \text{3d}_{yz} ,\, \text{3d}_{z^2} ,\, \text{3d}_{xz} ,\, \text{3d}_{x^2-y^2})$.
画出第一种表现形式的$\text{3d}$轨道, 也就是直接将量子数$(n, \ell, m) = (3, 2, m)$代入$\Psi_{n \ell m}$后得到的轨道, 其中磁量子数$m \in \mathbb{Z}$且$-2 \leq m \leq 2$.
根据以下变换关系, 画出第二种表现形式的$\text{3d}$轨道.
作为参考, 如下两张图分别是$\text{3d}_{-1}$轨道和$\text{3d}_{yz}$轨道:
画出第三章中其余三个$\text{sp}^3$杂化轨道的图像. 它们的形状应该与已有的轨道形状完全一致, 只是方向上有所区别.
$\text{2p}$轨道($\text{2p}_{x} ,\, \text{2p}_{y} ,\, \text{2p}_{z}$)是三个**等价**的原子轨道: 它们的形状完全一致, 只是朝向不同.
$\text{3d}$轨道($\text{3d}_{xy} ,\, \text{3d}_{yz} ,\, \text{3d}_{z^2} ,\, \text{3d}_{xz} ,\, \text{3d}_{x^2-y^2}$)大多是花瓣形的, 有四个分立的花瓣. 然而$\text{3d}_{z^2}$轨道并非如此, 它在$xOy$平面上有一个连通的环面, 让它显得格格不入.
于是就出现了$\text{3d}$轨道的另一种比较罕见的表现形式, 变换关系如下(注意轨道顺序):
\begin{equation} \begin{bmatrix} \Psi_{\text{3d, new1}} \\ \Psi_{\text{3d, new2}} \\ \Psi_{\text{3d, new3}} \\ \Psi_{\text{3d, new4}} \\ \Psi_{\text{3d, new5}} \end{bmatrix} = M \begin{bmatrix} \Psi_{\text{3d}_{z^2}} \\ \Psi_{\text{3d}_{xz}} \\ \Psi_{\text{3d}_{yz}} \\ \Psi_{\text{3d}_{x^2-y^2}} \\ \Psi_{\text{3d}_{xy}} \end{bmatrix} \end{equation}其中 \begin{equation} M = \begin{bmatrix} \sqrt{1 \over 5} & \sqrt{2 \over 5} & 0 & \sqrt{2 \over 5} & 0 \\ \sqrt{1 \over 5} & \sqrt{2 \over 5}\cos{2\pi \over 5} & \sqrt{2 \over 5}\sin{2\pi \over 5} & \sqrt{2 \over 5}\cos{4\pi \over 5} & \sqrt{2 \over 5}\sin{4\pi \over 5} \\ \sqrt{1 \over 5} & \sqrt{2 \over 5}\cos{4\pi \over 5} & \sqrt{2 \over 5}\sin{4\pi \over 5} & \sqrt{2 \over 5}\cos{8\pi \over 5} & \sqrt{2 \over 5}\sin{8\pi \over 5} \\ \sqrt{1 \over 5} & \sqrt{2 \over 5}\cos{6\pi \over 5} & \sqrt{2 \over 5}\sin{6\pi \over 5} & \sqrt{2 \over 5}\cos{12\pi \over 5} & \sqrt{2 \over 5}\sin{12\pi \over 5} \\ \sqrt{1 \over 5} & \sqrt{2 \over 5}\cos{8\pi \over 5} & \sqrt{2 \over 5}\sin{8\pi \over 5} & \sqrt{2 \over 5}\cos{16\pi \over 5} & \sqrt{2 \over 5}\sin{16\pi \over 5} \end{bmatrix} \end{equation}
(提示:下面是$\text{2p}$轨道的另一种表现形式)
\begin{equation} \begin{bmatrix} \Psi_{\text{2p, new1}} \\ \Psi_{\text{2p, new2}} \\ \Psi_{\text{2p, new3}} \end{bmatrix} = \begin{bmatrix} \sqrt{1 \over 3} & \sqrt{2 \over 3} & 0 \\ \sqrt{1 \over 3} & \sqrt{2 \over 3}\cos{2\pi \over 3} & \sqrt{2 \over 3}\sin{2\pi \over 3} \\ \sqrt{1 \over 3} & \sqrt{2 \over 3}\cos{4\pi \over 3} & \sqrt{2 \over 3}\sin{4\pi \over 3} \end{bmatrix} \begin{bmatrix} \Psi_{\text{2p}_{z}} \\ \Psi_{\text{2p}_{x}} \\ \Psi_{\text{2p}_{y}} \end{bmatrix} \end{equation}marching_cubes坐标互换的说明¶之前提到了marching_cubes的结果要做一次坐标互换, 似乎是scikit-image默认的坐标输出顺序有关.
验证是否需要互换/如何互换也比较简单, 用一个半主轴都不相等的椭球面即可. 例如下面这个椭球的半主轴长度分别为$8$, $2$和$5$.
\begin{equation} {x^2 \over 8^2} + {y^2 \over 2^2} + {z^2 \over 5^2} = 1 \end{equation}def ellipsoid(x, y, z):
return (x/8)**2 + (y/2)**2 + (z/5)**2
定义完之后, 就可以寻找等值面ellipsoid(x, y, z)=1
# 测试网格
test_limit = 10
test_n_points = 50
test_vec = np.linspace(-test_limit, test_limit, test_n_points)
Xt, Yt, Zt = np.meshgrid(test_vec, test_vec, test_vec)
test_fig, test_ax = new_fig_and_ax(plot_range=test_limit)
# 计算并处理椭球等值面(调整spacing & 平移至中心)
test_step = 2 * test_limit / (test_n_points - 1)
test_verts, test_faces, _, _ = marching_cubes(
ellipsoid(Xt, Yt, Zt),
level = 1,
spacing = (test_step, test_step, test_step),
)
test_verts -= test_limit
# 添加椭球等值面,带一定透明度以便于观察
test_mesh = Poly3DCollection(test_verts[test_faces], lw=0.1, alpha=0.5)
test_mesh.set_facecolor("pink")
test_mesh.set_edgecolor("grey")
test_ax.add_collection3d(test_mesh)
# 添加(8, 0, 0)、(0, 2, 0)、(0, 0, 5)作为参考
# 这些点都应该落在等值面(x/8)^2 + (y/2)^2 + (z/5)^2 = 1上
test_ax.scatter(8, 0, 0, color='red', s=20)
test_ax.scatter(0, 2, 0, color='green', s=20)
test_ax.scatter(0, 0, 5, color='blue', s=20)
plt.show()
三个点理应都落在上面, 然而只有$z$轴的$(0,0,5)$是正确的, 这是因为$xOy$平面上的两个轴方向颠倒了.
解决这个问题很简单,只需要互换test_verts的$x$坐标和$y$坐标即可.
# 互换所有顶点的x坐标和y坐标
test_verts[:,[0, 1]] = test_verts[:,[1, 0]]
# 重新绘图
new_test_fig, new_test_ax = new_fig_and_ax(plot_range=test_limit)
test_mesh = Poly3DCollection(test_verts[test_faces], lw=0.1, alpha=0.5)
test_mesh.set_facecolor("yellow")
test_mesh.set_edgecolor("grey")
new_test_ax.add_collection3d(test_mesh)
new_test_ax.scatter(8, 0, 0, color='red', s=20)
new_test_ax.scatter(0, 2, 0, color='green', s=20)
new_test_ax.scatter(0, 0, 5, color='blue', s=20)
# 可以发现这次三个点都落在等值面上,问题解决
plt.show()